---
title: "Analysis"
author: "Ksiazkiewicz et al."
date: "2024-02-27"
output:
  word_document: default
  html_document: default
  pdf_document: default
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = FALSE)
library(readr)
library(tidyverse)
library(plyr)
library(psych)
library(ggplot2)
library(texreg)


data2=read.csv("~/[adjust the path as needed]/df2022.csv")
data=read.csv("~/[adjust the path as needed]/df2023.csv")
```

Study 2 Factor Analysis (Genetic Determinism, Risk, Concerned, Trust, Experience, Willingness)
Study 3 Factor Analysis (Risk, Concerned, Trust, Experience, Willingness, GSPQ)
Study 3 Item Response Theory


## Factor Analysis: Genetic Determinism (Study 2)
```{r}
#genetic determinism
gendeter<-data2[,c("Q56","Q57", "Q61", "Q62", "Q63", "Q58","Q59", "Q60")]
gendeter<-na.omit(gendeter)
pc.gendeter<-princomp(gendeter, cor=TRUE)

#eigenvectors
pc.gendeter$loadings
#eigenvalues
pc.gendeter$sdev *pc.gendeter$sdev

round(cor(gendeter, pc.gendeter$scores), 2)
screeplot(pc.gendeter, type = "l", main="screeplot") 
abline(1,0, col= 'red', lty = 2)

data2$gendeter=(data2$Q56+data2$Q57+data2$Q58+data2$Q59+data2$Q60+data2$Q61+data2$Q62+data2$Q63)/8

ps3<-principal(gendeter, nfactors = 2, rotate = "none")
ps3 

#PC1
#Q56: PC1 ●	I think the chief reason why parents and children are so alike in behavior and character is that they possess a shared genetic inheritance.
#Q57: PC1 ●	I believe that children inherit many of their personal traits from their parents
#Q59: PC1 ●	I think that the upbringing by parents and the social environment have far greater significance for the development of abilities and personal traits than genetic predispositions. (R)
#Q61: PC1 ●	I am convinced that the analysis of the genetic predispositions of an embryo allows good predictions as to which characteristic and abilities the child will develop.
#Q62: PC1 ●	I think that twins, because of the identical genetic predispositions, will be very similar in their behavior even if they were adopted and raised in different families
#Q63: PC1 ●	I believe that an analysis of my genetic predispositions will allow a trained scientist to predict many of my abilities and traits without having any personal knowledge of me.
data2$gendeter1=(data2$Q56+data2$Q57+data2$Q59+data2$Q61+data2$Q62+data2$Q63)/6
#PC2
#Q58: PC2 ●	I am convinced that very few behavioral traits of humans can be traced back to their genes. (R)
#Q60: PC2 ●	I think that genetic predispositions have little influence on a person’s personality characteristics. (R)
data2$gendeter2=(data2$Q58+data2$Q60)/2
```


## Factor Analysis: Risk (Study 2)
```{r}
#risk (higher is risk-taking)
risk1<-data2[,c("Q76","Q77", "Q78", "Q79", "Q80", "Q81","Q82", "Q83")]
risk1<-na.omit(risk1)
pc.risk1<-princomp(risk1, cor=TRUE)

#eigenvectors
pc.risk1$loadings
#eigenvalues
pc.risk1$sdev *pc.risk1$sdev

round(cor(risk1, pc.risk1$scores), 2)
screeplot(pc.risk1, type = "l", main="screeplot") 
abline(1,0, col= 'red', lty = 2)

data2$risk=(data2$Q76+data2$Q77+data2$Q78+data2$Q79+data2$Q80+data2$Q81+data2$Q82+data2$Q83)/8
```
## Factor Analysis: Risk (Study 3)
```{r}
#risk (higher is risk-taking)
risk2<-data[,c("Q76","Q77", "Q78", "Q79", "Q80", "Q81","Q82", "Q83")]
risk2<-na.omit(risk2)
pc.risk2<-princomp(risk2, cor=TRUE)

#eigenvectors
pc.risk2$loadings
#eigenvalues
pc.risk2$sdev *pc.risk2$sdev

round(cor(risk2, pc.risk2$scores), 2)
screeplot(pc.risk2, type = "l", main="screeplot") 
abline(1,0, col= 'red', lty = 2)

data$risk=(data$Q76+data$Q77+data$Q78+data$Q79+data$Q80+data$Q81+data$Q82+data$Q83)/8
```


## Factor Analysis: Concerned (Study 2)
```{r}
#concerned
concern1<-data2[,c("Q17_1", "Q17_2", "Q17_3", "Q18_1", "Q18_2", "Q18_3", "Q18_4", "Q18_5", "Q19_1", "Q19_2", "Q19_3", "Q19_4", "Q19_5", "Q19_6")]
concern1<-na.omit(concern1)
pc.concern1<-princomp(concern1, cor=TRUE)

#eigenvectors
pc.concern1$loadings
#eigenvalues
pc.concern1$sdev *pc.concern1$sdev

round(cor(concern1, pc.concern1$scores), 2)
screeplot(pc.concern1, type = "l", main="screeplot") 
abline(1,0, col= 'red', lty = 2)

data2$concern=data2$Q17_1+data2$Q17_2+data2$Q17_3+data2$Q18_1+data2$Q18_2+data2$Q18_3+data2$Q18_4+data2$Q18_5+data2$Q19_1+data2$Q19_2+data2$Q19_3+data2$Q19_4+data2$Q19_5+data2$Q19_6
```
## Factor Analysis: Concerned (Study 3)
```{r}
#concerned
concern2<-data[,c("Q17_1", "Q17_2", "Q17_3", "Q18_1", "Q18_2", "Q18_3", "Q18_4", "Q18_5", "Q19_1", "Q19_2", "Q19_3", "Q19_4", "Q19_5", "Q19_6")]
concern2<-na.omit(concern2)
pc.concern2<-princomp(concern2, cor=TRUE)
#eigenvectors
pc.concern2$loadings
#eigenvalues
pc.concern2$sdev *pc.concern2$sdev
round(cor(concern2, pc.concern2$scores), 2)
screeplot(pc.concern2, type = "l", main="screeplot")
abline(1,0, col= 'red', lty = 2)

data$concern=(data$Q17_1+data$Q17_2+data$Q17_3+data$Q18_1+data$Q18_2+data$Q18_3+data$Q18_4+data$Q18_5+data$Q19_1+data$Q19_2+data$Q19_3+data$Q19_4+data$Q19_5+data$Q19_6)/14
```


## Factor Analysis: Trust (Study 2)
```{r}
#trust
trust1<-data2[,c("Q84_1", "Q84_2", "Q84_3", "Q84_7", "Q84_5", "Q84_6")]
trust1<-na.omit(trust1)
pc.trust1<-princomp(trust1, cor=TRUE)

#eigenvectors
pc.trust1$loadings
#eigenvalues
pc.trust1$sdev *pc.trust1$sdev

round(cor(trust1, pc.trust1$scores), 2)

screeplot(pc.trust1, type = "l", main = "Screeplot")
abline(h = 1, col = 'red', lty = 2)

data2$trust=(data2$Q84_1+data2$Q84_2+data2$Q84_3+data2$Q84_7+data2$Q84_5+data2$Q84_6)/6
data2$trust3=(data2$Q84_1+data2$Q84_2+data2$Q84_3)/3
```
## Factor Analysis: Trust (Study 3)
```{r}
#trust
trust2<-data[,c("Q84_1", "Q84_2", "Q84_3", "Q84_7", "Q84_5", "Q84_6", "Q84_8", "Q84_9")]
trust2<-na.omit(trust2)
pc.trust2<-princomp(trust2, cor=TRUE)

#eigenvectors
pc.trust2$loadings
#eigenvalues
pc.trust2$sdev *pc.trust2$sdev

round(cor(trust2, pc.trust2$scores), 2)
screeplot(pc.trust2, type = "l", main="screeplot") 
abline(1,0, col= 'red', lty = 2)

data$trust=(data$Q84_1+data$Q84_2+data$Q84_3+data$Q84_7+data$Q84_5+data$Q84_6+data$Q84_8+data$Q84_9)/8
```


## Factor Analysis: Experience (Study 2)
```{r}
#Experience
experience1<-data2[,c("exp1", "exp2","exp3")]
experience1<-na.omit(experience1)
pc.experience1<-princomp(experience1, cor=TRUE)
#eigenvectors
pc.experience1$loadings
#eigenvalues
pc.experience1$sdev *pc.experience1$sdev
round(cor(experience1, pc.experience1$scores), 2)
screeplot(pc.experience1, type = "l", main="screeplot") 
abline(1,0, col= 'red', lty = 2)

data2$exp=data2$exp1+data2$exp2+data2$exp3
data2$exp=as.factor(data2$exp)
```
## Factor Analysis: Experience (Study 3)
```{r}
#experience
experience2<-data[,c("exp1", "exp2","exp3")]
experience2<-na.omit(experience2)
pc.experience2<-princomp(experience2, cor=TRUE)
#eigenvectors
pc.experience2$loadings
#eigenvalues
pc.experience2$sdev *pc.experience2$sdev
round(cor(experience2, pc.experience2$scores), 2)
screeplot(pc.experience2, type = "l", main="screeplot") 
abline(1,0, col= 'red', lty = 2)

data$exp=data$exp1+data$exp2+data$exp3
data$exp=as.factor(data$exp)
```


## Factor Analysis: Willingness (Study 2)
```{r}
#willingness
willingness1<-data2[,c("will2", "will3","will4")]
willingness1<-na.omit(willingness1)
pc.willingness1<-princomp(willingness1, cor=TRUE)
#eigenvectors
pc.willingness1$loadings
#eigenvalues
pc.willingness1$sdev *pc.willingness1$sdev
round(cor(willingness1, pc.willingness1$scores), 2)
screeplot(pc.willingness1, type = "l", main="screeplot") 
abline(1,0, col= 'red', lty = 2)
data2$will=(data2$will2+data2$will3+data2$will4)/3
```
## Factor Analysis: Willingness (Study 3)
```{r}
#willingness
willingness2<-data[,c("will1", "will2", "will3","will4")]
willingness2<-na.omit(willingness2)
pc.willingness2<-princomp(willingness2, cor=TRUE)
#eigenvectors
pc.willingness2$loadings
#eigenvalues
pc.willingness2$sdev *pc.willingness2$sdev
round(cor(willingness2, pc.willingness2$scores), 2)
screeplot(pc.willingness2, type = "l", main="screeplot") 
abline(1,0, col= 'red', lty = 2)
data$will=(data$will1+data$will2+data$will3+data$will4)/4
```



```{r}
#GSPQ
#factor analysis
gspq<-data[,c("GSPQ_1", "GSPQ_13", "GSPQ_20", "GSPQ_28", "GSPQ_36", "GSPQ_2", "GSPQ_14", "GSPQ_29", "GSPQ_37", "GSPQ_3", "GSPQ_15", "GSPQ_21", "GSPQ_30", "GSPQ_38", "GSPQ_4", "GSPQ_16", "GSPQ_22", "GSPQ_39", "GSPQ_5", "GSPQ_32", "GSPQ_50", "GSPQ_53", "GSPQ_54", "GSPQ_55", "GSPQ_56", "GSPQ_58", "GSPQ_64", "GSPQ_69", "GSPQ_6", "GSPQ_17", "GSPQ_24", "GSPQ_33", "GSPQ_40", "GSPQ_51", "GSPQ_7", "GSPQ_18", "GSPQ_25", "GSPQ_34", "GSPQ_8", "GSPQ_19", "GSPQ_27", "GSPQ_35", "GSPQ_42", "GSPQ_52", "GSPQ_9", "GSPQ_60", "GSPQ_10", "GSPQ_11", "GSPQ_12", "GSPQ_59", "GSPQ_23", "GSPQ_26", "GSPQ_31", "GSPQ_41", "GSPQ_43", "GSPQ_44", "GSPQ_45", "GSPQ_46", "GSPQ_47", "GSPQ_48", "GSPQ_49", "GSPQ_57", "GSPQ_61", "GSPQ_62", "GSPQ_63", "GSPQ_65", "GSPQ_66", "GSPQ_67", "GSPQ_68", "GSPQ_70", "GSPQ_71", "GSPQ_72", "GSPQ_73")]
gspq<-na.omit(gspq)
pc.gspq<-princomp(gspq, cor=TRUE)

#eigenvectors
pc.gspq$loadings
#eigenvalues
pc.gspq$sdev *pc.gspq$sdev

round(cor(gspq, pc.gspq$scores), 2)
screeplot(pc.gspq, type = "l", main="screeplot") 
abline(1,0, col= 'red', lty = 2)

ps_gspq<-principal(gspq, nfactors = 3, rotate = "none")
ps_gspq

round(cor(gspq, ps_gspq$scores), 2)


data$pc1=(data$GSPQ_1+data$GSPQ_13+data$GSPQ_20+data$GSPQ_28+data$GSPQ_36+data$GSPQ_2+data$GSPQ_14+data$GSPQ_29+data$GSPQ_37+data$GSPQ_6+data$GSPQ_17+data$GSPQ_24+data$GSPQ_33+data$GSPQ_40+data$GSPQ_51+data$GSPQ_7+data$GSPQ_18+data$GSPQ_25+data$GSPQ_34+data$GSPQ_8+data$GSPQ_19+data$GSPQ_27+data$GSPQ_35+data$GSPQ_42+data$GSPQ_52+data$GSPQ_9+data$GSPQ_60+data$GSPQ_11+data$GSPQ_65+data$GSPQ_66)/30

data$pc2=(data$GSPQ_3+data$GSPQ_15+data$GSPQ_21+data$GSPQ_30+data$GSPQ_38+data$GSPQ_4+data$GSPQ_16+data$GSPQ_22+data$GSPQ_39+data$GSPQ_50+data$GSPQ_55+data$GSPQ_56+data$GSPQ_58+data$GSPQ_64+data$GSPQ_69+data$GSPQ_23+data$GSPQ_26+data$GSPQ_43+data$GSPQ_45+data$GSPQ_46+data$GSPQ_47+data$GSPQ_57+data$GSPQ_61+data$GSPQ_67+data$GSPQ_68+data$GSPQ_70+data$GSPQ_71+data$GSPQ_72+data$GSPQ_73)/29

data$pc3=(data$GSPQ_5+data$GSPQ_32+data$GSPQ_53+data$GSPQ_54+data$GSPQ_31+data$GSPQ_41+data$GSPQ_44+data$GSPQ_48+data$GSPQ_49+data$GSPQ_63+(6-data$GSPQ_62)+(6-data$GSPQ_10)+(6-data$GSPQ_12)+(6-data$GSPQ_59))/14
```


```{r}
pc1=data$pc1
pc2=data$pc2
pc3=data$pc3
cor.test(pc1, pc2)#-0.09
ggplot(data, aes(x=pc1, y = pc2)) + geom_point() + geom_smooth(method=lm, se=FALSE) + geom_text(data=data.frame(x=c(1.5,1.5), y=c(5,5)), aes(x=x, y=y, label="Cor=-0.09", size=10), inherit.aes = FALSE, show.legend = F)
cor.test(pc2, pc3)#0.21
ggplot(data, aes(x=pc2, y = pc3)) + geom_point() + geom_smooth(method=lm, se=FALSE) + geom_text(data=data.frame(x=c(1.5,1.5), y=c(5,5)), aes(x=x, y=y, label="Cor=0.21", size=10), inherit.aes = FALSE, show.legend = F)
cor.test(pc1, pc3) #-0.14
ggplot(data, aes(x=pc1, y = pc3)) + geom_point() + geom_smooth(method=lm, se=FALSE) + geom_text(data=data.frame(x=c(1.5,1.5), y=c(5,5)), aes(x=x, y=y, label="Cor=-0.14", size=10), inherit.aes = FALSE, show.legend = F)
```


```{r}
#configural invariance
library(lavaan)
cfa=gspq<-data[,c("GSPQ_1", "GSPQ_13", "GSPQ_20", "GSPQ_28", "GSPQ_36", "GSPQ_2", "GSPQ_14", "GSPQ_29", "GSPQ_37", "GSPQ_3", "GSPQ_15", "GSPQ_21", "GSPQ_30", "GSPQ_38", "GSPQ_4", "GSPQ_16", "GSPQ_22", "GSPQ_39", "GSPQ_5", "GSPQ_32", "GSPQ_50", "GSPQ_53", "GSPQ_54", "GSPQ_55", "GSPQ_56", "GSPQ_58", "GSPQ_64", "GSPQ_69", "GSPQ_6", "GSPQ_17", "GSPQ_24", "GSPQ_33", "GSPQ_40", "GSPQ_51", "GSPQ_7", "GSPQ_18", "GSPQ_25", "GSPQ_34", "GSPQ_8", "GSPQ_19", "GSPQ_27", "GSPQ_35", "GSPQ_42", "GSPQ_52", "GSPQ_9", "GSPQ_60", "GSPQ_10", "GSPQ_11", "GSPQ_12", "GSPQ_59", "GSPQ_23", "GSPQ_26", "GSPQ_31", "GSPQ_41", "GSPQ_43", "GSPQ_44", "GSPQ_45", "GSPQ_46", "GSPQ_47", "GSPQ_48", "GSPQ_49", "GSPQ_57", "GSPQ_61", "GSPQ_62", "GSPQ_63", "GSPQ_65", "GSPQ_66", "GSPQ_67", "GSPQ_68", "GSPQ_70", "GSPQ_71", "GSPQ_72", "GSPQ_73", "Race")]
cfa<-na.omit(cfa)

cfa$Race = as.factor(cfa$Race)

cfa <- cfa %>%
  mutate(
    GSPQ_10 = 6 - GSPQ_10,
    GSPQ_12 = 6 - GSPQ_12,
    GSPQ_59 = 6 - GSPQ_59,
    GSPQ_62 = 6 - GSPQ_62
  )

model <- '
  # First factor
  pc1 =~ GSPQ_1 + GSPQ_13 + GSPQ_20 + GSPQ_28 + GSPQ_36 + GSPQ_2 + GSPQ_14 + GSPQ_29 + GSPQ_37 +
         GSPQ_6 + GSPQ_17 + GSPQ_24 + GSPQ_33 + GSPQ_40 + GSPQ_51 + GSPQ_7 + GSPQ_18 + GSPQ_25 + GSPQ_34 +
         GSPQ_8 + GSPQ_19 + GSPQ_27 + GSPQ_35 + GSPQ_42 + GSPQ_52 + GSPQ_9 + GSPQ_60 + GSPQ_11 + GSPQ_65 + GSPQ_66

  # Second factor
  pc2 =~ GSPQ_3 + GSPQ_15 + GSPQ_21 + GSPQ_30 + GSPQ_38 + GSPQ_4 + GSPQ_16 + GSPQ_22 + GSPQ_39 +
         GSPQ_50 + GSPQ_55 + GSPQ_56 + GSPQ_58 + GSPQ_64 + GSPQ_69 + GSPQ_23 + GSPQ_26 + GSPQ_43 + GSPQ_45 +
         GSPQ_46 + GSPQ_47 + GSPQ_57 + GSPQ_61 + GSPQ_67 + GSPQ_68 + GSPQ_70 + GSPQ_71 + GSPQ_72 + GSPQ_73

  # Third factor (with reverse coded items)
  pc3 =~ GSPQ_5 + GSPQ_32 + GSPQ_53 + GSPQ_54 + GSPQ_31 + GSPQ_41 + GSPQ_44 + GSPQ_48 + GSPQ_49 + GSPQ_63 +
         GSPQ_10 + GSPQ_12 + GSPQ_59 + GSPQ_62
'

# Fit model by group
results <- list()
races <- levels(cfa$Race)

for(race in races) {
  group_data <- subset(cfa, Race == race)
  fit <- cfa(model, data=group_data, missing="ml")
  results[[race]] <- summary(fit, fit.measures=TRUE)
  print(paste("Results for Race:", race))
  print(fitMeasures(fit, c("cfi", "tli", "rmsea", "srmr")))}
```





```{r}
gspq1<-data[,c("GSPQ_1", "GSPQ_13", "GSPQ_20", "GSPQ_28", "GSPQ_36", "GSPQ_2", "GSPQ_14", "GSPQ_29", "GSPQ_37", "GSPQ_6", "GSPQ_17", "GSPQ_24", "GSPQ_33", "GSPQ_40", "GSPQ_51", "GSPQ_7", "GSPQ_18", "GSPQ_25", "GSPQ_34", "GSPQ_8", "GSPQ_19", "GSPQ_27", "GSPQ_35", "GSPQ_42", "GSPQ_52", "GSPQ_9", "GSPQ_60", "GSPQ_11", "GSPQ_65", "GSPQ_66")]
     
gspq2<-data[,c("GSPQ_3", "GSPQ_15", "GSPQ_21", "GSPQ_30", "GSPQ_38", "GSPQ_4", "GSPQ_16", "GSPQ_22", "GSPQ_39", "GSPQ_50", "GSPQ_55", "GSPQ_56", "GSPQ_58", "GSPQ_64", "GSPQ_69", "GSPQ_23", "GSPQ_26", "GSPQ_43", "GSPQ_45", "GSPQ_46", "GSPQ_47", "GSPQ_57", "GSPQ_61", "GSPQ_67", "GSPQ_68", "GSPQ_70", "GSPQ_71", "GSPQ_72", "GSPQ_73")]

gspq3<-data[,c("GSPQ_5", "GSPQ_32",  "GSPQ_53", "GSPQ_54", "GSPQ_31", "GSPQ_41",  "GSPQ_44",  "GSPQ_48", "GSPQ_49", "GSPQ_63", "GSPQ_10", "GSPQ_12", "GSPQ_59", "GSPQ_62")] 
gspq3= gspq3 %>% mutate_at(c("GSPQ_10", "GSPQ_12", "GSPQ_59", "GSPQ_62"), funs(recode(., "5"=1, "4"=2, "3"=3, "2"=4, "1"=5)))

alpha(gspq1)$total #0.94
alpha(gspq2)$total #0.92
alpha(gspq3)$total #0.81

gspq1<-na.omit(gspq1)
gspq2<-na.omit(gspq2)
gspq3<-na.omit(gspq3)
```
# Model fit
```{r}
library(mirt)
mod1 <- (mirt(gspq1, 1, verbose = FALSE, itemtype = 'graded', SE = TRUE))
M2(mod1, type = "C2", calcNULL = FALSE)
mod2 <- (mirt(gspq2, 1, verbose = FALSE, itemtype = 'graded', SE = TRUE))
M2(mod2, type = "C2", calcNULL = FALSE)
mod3 <- (mirt(gspq3, 1, verbose = FALSE, itemtype = 'graded', SE = TRUE))
M2(mod3, type = "C2", calcNULL = FALSE)
```

# Item fit
```{r}
itemfit(mod1)
itemfit(mod2)
itemfit(mod3)
```

# IRT parameters:
```{r}
library(ltm)
grm(gspq1)
grm(gspq2)
grm(gspq3)

IRT_parms1 <- coef(mod1, IRTpars = TRUE, simplify = TRUE)
IRT_parms1$items
IRT_parms2 <- coef(mod2, IRTpars = TRUE, simplify = TRUE)
IRT_parms2$items
IRT_parms3 <- coef(mod3, IRTpars = TRUE, simplify = TRUE)
IRT_parms3$items

# Partial Credit Model
results.gpcm <- mirt(data=gspq1, itemtype="gpcm", SE=TRUE, verbose=FALSE)
coef.gpcm <- coef(results.gpcm, IRTpars=TRUE, simplify=TRUE)
items.gpcm <- as.data.frame(coef.gpcm$items)
print(items.gpcm)

results.gpcm <- mirt(data=gspq2, itemtype="gpcm", SE=TRUE, verbose=FALSE)
coef.gpcm <- coef(results.gpcm, IRTpars=TRUE, simplify=TRUE)
items.gpcm <- as.data.frame(coef.gpcm$items)
print(items.gpcm)

results.gpcm <- mirt(data=gspq3, itemtype="gpcm", SE=TRUE, verbose=FALSE)
coef.gpcm <- coef(results.gpcm, IRTpars=TRUE, simplify=TRUE)
items.gpcm <- as.data.frame(coef.gpcm$items)
print(items.gpcm) #44 vs 63

# Rating Scale Model
results.rsm <- mirt(data=gspq1, itemtype="rsm", verbose=FALSE)
coef.rsm <- coef(results.rsm, simplify=TRUE)
items.rsm <- as.data.frame(coef.rsm$items)
print(items.rsm)
```

# IRT Plots
```{r}
# Item information curves
plot(mod1, type='infotrace', which.item = c(1:30), facet_items=T, as.table = TRUE, auto.key=list(points=F, lines=T, columns=1, space = 'right', cex = .8), main="")

plot(mod2, type='infotrace', which.item = c(1:29), facet_items=T, as.table = TRUE, auto.key=list(points=F, lines=T, columns=1, space = 'right', cex = .8), main="")

plot(mod3, type='infotrace', which.item = c(1:14), facet_items=T, as.table = TRUE, auto.key=list(points=F, lines=T, columns=1, space = 'right', cex = .8), theta_lim = c(-3, 3), main="")
```

#GSPQ Final 16 Study 3
```{r}
data$pc1_final=(data$GSPQ_7+data$GSPQ_18+data$GSPQ_27+data$GSPQ_29+data$GSPQ_42+data$GSPQ_65)/6
data$pc2_final=(data$GSPQ_3+data$GSPQ_16+data$GSPQ_38+data$GSPQ_39+data$GSPQ_67)/5
data$pc3_final=(data$GSPQ_5+data$GSPQ_32+data$GSPQ_49+data$GSPQ_53+data$GSPQ_54)/5
```
